# 用于对训练数据和流动电位仪做数据过滤，并输出过滤后的结果。
from typing import Literal
from numpy.lib.stride_tricks import sliding_window_view
import pandas as pd
import numpy as np
from filterpy.kalman import KalmanFilter

# 封装滤波算法的类
class Filter:
    def __init__(self, data, window_size=10, threshold=0.1):
        self.data = data
        self.window_size = window_size
        self.threshold = threshold

    def moving_average(self):
        """
        移动平均滤波
        """
        result = self.data.rolling(window=self.window_size).mean()
        return result
    
    def adaptive_moving_average(self,data, window_size, threshold=0.5):
    # 自适应滑动平均
        filtered = np.copy(data)
        window = window_size
        for i in range(len(data)):
            if i >= window:
                window = window_size + int(np.abs(data[i] - np.mean(data[i-window:i])) / threshold)
            filtered[i] = np.mean(data[max(0, i-window):i+1])
        return filtered

    # 中值滤波
    def median_filter(self, data, window_size):
            filtered = np.zeros_like(data)
            half_window = window_size // 2 # 窗口大小的一半 //表示除法，返回商的整数部分
            for i in range(len(data)):
                window = data[max(0, i-half_window):min(len(data), i+half_window+1)]
                filtered[i] = np.median(window)
            return filtered
  


    def threshold_filter(self):
        """
        阈值滤波
        """
        result = self.data.copy()
        result[result < self.threshold] = np.nan
        return result

    # 均值滤波
    def mean_filter(self):
        """
        均值滤波
        """
        result = self.data.rolling(window=self.window_size).mean()
        return result
    # 低通滤波器
    def low_pass_filter(self, cutoff=0.1):
        """
        低通滤波器
        """
        result = self.data.copy()
        result[result > cutoff] = cutoff
        return result
    # 高通滤波器
    def high_pass_filter(self, cutoff=0.1):
        """
        高通滤波器
        """
        result = self.data.copy()
        result[result < cutoff] = cutoff
        return result
    # 双边滤波器
    def band_pass_filter(self, low_cutoff=0.1, high_cutoff=0.5):
        """
        双边滤波器
        """
        result = self.data.copy()
        result[result < low_cutoff] = np.nan
        result[result > high_cutoff] = np.nan
        return result
    # 双边带通滤波器
    def band_stop_filter(self, low_cutoff=0.1, high_cutoff=0.5):
        """
        双边带通滤波器
        """
        result = self.data.copy()
        result[result > low_cutoff] = np.nan
        result[result < high_cutoff] = np.nan
        return result
    

    def particle_filter(
        observations,                # 观测数据（形状 [num_timesteps, obs_dim]）
        transition_func,             # 状态转移函数（输入粒子数组，输出下一状态粒子数组）
        observation_func,            # 观测函数（输入粒子数组，输出预测观测数组）
        initial_particles,           # 初始粒子数组（形状 [num_particles, state_dim]）
        process_noise_cov,           # 过程噪声协方差矩阵（形状 [state_dim, state_dim]）
        observation_noise_cov,       # 观测噪声协方差矩阵（形状 [obs_dim, obs_dim]）
        num_particles,               # 粒子数量（整数）
        resample_threshold=0.5       # 重采样阈值（ESS < 粒子数 * threshold 时触发）
    ):
        """
        执行粒子滤波算法，返回状态估计结果
        
        参数:
            observations: 观测数据数组（时间步 × 观测维度）
            transition_func: 状态转移函数（粒子数组 → 下一状态粒子数组）
            observation_func: 观测函数（粒子数组 → 预测观测数组）
            initial_particles: 初始粒子数组（粒子数 × 状态维度）
            process_noise_cov: 过程噪声协方差矩阵
            observation_noise_cov: 观测噪声协方差矩阵
            num_particles: 粒子数量
            resample_threshold: 重采样触发阈值（默认0.5，即ESS < 粒子数×0.5时重采样）
        
        返回:
            estimates: 状态估计数组（时间步 × 状态维度）
        """
        num_timesteps = observations.shape[0]
        state_dim = initial_particles.shape[1]
        obs_dim = observations.shape[1]
        estimates = np.zeros((num_timesteps, state_dim))  # 存储状态估计
        
        # 初始化粒子、权重、ESS
        particles = initial_particles.copy()
        weights = np.ones(num_particles) / num_particles
        ess = num_particles  # 初始有效样本大小
        
        for t in range(num_timesteps):
            # ---------------------- 预测步骤 ----------------------
            # 传播粒子并添加过程噪声
            particles = transition_func(particles)  # 应用状态转移模型
            process_noise = np.random.multivariate_normal(
                mean=np.zeros(state_dim),
                cov=process_noise_cov,
                size=num_particles
            )
            particles += process_noise  # 添加过程噪声
            
            # ---------------------- 更新步骤 ----------------------
            # 计算预测观测值（当前粒子对应的观测）
            predicted_obs = observation_func(particles)  # 形状 [num_particles, obs_dim]
            
            # 计算每个粒子的观测似然（高斯分布）
            likelihoods = np.zeros(num_particles)
            for i in range(num_particles):
                residual = observations[t] - predicted_obs[i]
                d = obs_dim
                det_R = np.linalg.det(observation_noise_cov)
                inv_R = np.linalg.inv(observation_noise_cov)
                exponent = -0.5 * residual.T @ inv_R @ residual
                likelihood = (1 / (np.power(2 * np.pi, d/2) * np.sqrt(det_R))) * np.exp(exponent)
                likelihoods[i] = likelihood
            
            # 更新权重（乘以似然并归一化）
            weights *= likelihoods
            weights += 1e-300  # 防止除以0
            weights /= np.sum(weights)
            
            # 计算有效样本大小（ESS）
            ess = 1.0 / np.sum(weights ** 2)
            
            # ---------------------- 重采样步骤 ----------------------
            if ess < num_particles * resample_threshold:
                # 系统重采样
                cumulative_weights = np.cumsum(weights)
                u = np.random.uniform(0, 1/num_particles)
                indices = np.zeros(num_particles, dtype=int)
                j = 0
                for i in range(num_particles):
                    while cumulative_weights[j] < u:
                        j += 1
                    indices[i] = j
                    u += 1/num_particles
                particles = particles[indices]  # 替换为高权重粒子
                weights = np.ones(num_particles) / num_particles  # 重置权重
                ess = num_particles  # 重采样后ESS=粒子数
            
            # ---------------------- 状态估计 ----------------------
            estimates[t] = np.mean(particles, axis=0)  # 粒子均值作为状态估计
        
        return estimates
    # 卡尔曼滤波器
    def kalman_filter(self, data, process_noise=0.1, measurement_noise=1.0):
        """
        执行卡尔曼滤波算法，返回状态估计结果
        
        参数:
            observations: 观测数据数组（时间步 × 观测维度）
            initial_state: 初始状态向量（状态维度）
            state_transition: 状态转移矩阵 F（状态维度 × 状态维度）
            observation_matrix: 观测矩阵 H（观测维度 × 状态维度）
            process_noise_cov: 过程噪声协方差 Q（状态维度 × 状态维度）
            observation_noise_cov: 观测噪声协方差 R（观测维度 × 观测维度）
        
        返回:
            estimates: 状态估计数组（时间步 × 状态维度）
        """
        kf = KalmanFilter(dim_x=1, dim_z=1)
        kf.x = np.array([data[0]])  # 初始状态
        kf.F = np.array([[1.]])      # 状态转移矩阵
        kf.H = np.array([[1.]])      # 观测矩阵
        kf.P *= 1000.                # 初始协方差
        kf.R = measurement_noise     # 观测噪声
        kf.Q = process_noise         # 过程噪声

        filtered = []
        for z in data:
            kf.predict()
            kf.update(z)
            filtered.append(kf.x[0])
        return np.array(filtered)
    # def Ekalman_filer()
        





